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Abstract 

We propose a method for pricing American options whose pay-off depends on the moving 
average of the underlying asset price. The method uses a finite dimensional approximation of 
the infinite-dimensional dynamics of the moving average process based on a truncated Laguerre 
series expansion. The resulting problem is a finite-dimensional optimal stopping problem, which 
we propose to solve with a least squares Monte Carlo approach. We analyze the theoretical 
convergence rate of our method and present numerical results in the Black-Scholes framework. 
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1 Introduction 

We introduce a new method to value American options whose payoff at exercise depends on the 
moving average of the underlying asset price. The simplest example (sometimes known as surge 
option) is a variable strike call or put, whose strike is adjusted daily to the moving average of 
the underlying asset over a certain fixed-length period preceding the current date. American-style 
options on moving average are widely used in energy markets. In gas markets, for example, these 
options are known as indexed Swing options and allow the holder to purchase an amount of gas at 
a strike price, which is indexed on moving averages of various oil-prices: typically gas oil and fuel 
oil prices are averaged over the last 6 months and delayed in time with a 1 month lag. 

We shall denote by X the moving average of the underlying S over a time window with fixed 
length 5 > 0: 

i r* 

X t = - S u du, Vt > 5. (1) 

o Jt-S 
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The process X follows the dynamics 



dX t 



\ (S t - Sts) dt, Vt > S. 



This shows in particular that even if S is Markovian, the process (S, X) is not: it is, in general, 
impossible for any finite n to find n processes X , . . . , X n such that (S, X, X , . . . , X n ) are jointly 
Markovian. This property makes the pricing of the moving window options with early exercise 
a challenging problem both from the theoretical and the numerical viewpoint. In a continuous 
time framework the problem is infinite dimensional, and in a discrete time framework (pricing of a 
Bermudan option instead of an American option) there is a computational challenge, due to high 
dimensionality: the dimension is equal to the number of time steps within the averaging window. 
This in particular makes it difficult to compute the conditional expectations involved in the optimal 
exercise rule. 

The problem of pricing moving average American options should not be confused with a much 
simpler problem of pricing Asian American options with a fixed start averaging window, where the 
payoff depends on 



It is well-known (see for example Wilmott and al. (25) ) that in this case, adding a dimension to the 
problem allows to derive a finite dimensional Markovian formulation. On the other hand, partial 
average Asian options of European style can be easily valued (see for example Shreeve [22]). If the 
averaging period has a length S > 0, then on [T — 5, T] the option value is given by the price of 
the corresponding Asian option and on [0, T — 8] it solves a European style PDE with appropriate 
terminal and boundary conditions. 

In this paper, we propose a method for pricing moving average American options based on 
a finite dimensional approximation of the infinite-dimensional dynamics of the moving average 
process. The approximation is based on a truncated expansion of the weighting measure used for 
averaging in a series involving Laguerre polynomials. This technique has long been used in signal 
processing (see for example [T7]) but is less known in the context of approximation of stochastic 
systems. We compute the rate of convergence of our method as function of the number of terms in 
the series. The resulting problem is then a finite-dimensional optimal stopping problem, which we 
propose to solve with a Monte Carlo Longstaff and Schwartz-type approach. Numerical results are 
presented for moving average options in the Black-Scholcs framework. 

In the literature, very few articles discuss moving average options with early exercise feature 
[U |U El [TTJ [T2"]. A common approach (see e.g., Broadie and Cao [I]) is to use the least squares 
Monte Carlo, computing the conditional expectation estimators through regressions on polynomial 
functions of the current values of the underlying price and its moving average. Since the future 
evolution of the moving average depends on the entire history of the price process between t — 5 
and t, this approach introduces a bias. In our numerical examples we compare this approach to our 
results, and find that for standard moving average American options the error is not so large (less 
than 1% for the examples we took), which justifies the use of this approach for practical purposes 
in spite of its suboptimality. For moving average American options with time delay, whose payoff 
depends on the average of the price between dates t — Si and t — <5 2 , < S 2 < Si , the suboptimal 
approximation leads to a bias of up to 11% of the option's price in our examples. 

Bilger [3] uses a regression based approach in the discrete-time setting to compute the con- 
ditional expectations considering that the state vector is composed of the underlying price, its 
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moving average and additional partial averages of the price over the rolling period. Their number 
is computed heuristically and as it tends to the number of time steps within the rolling period, the 
computed price tends to the true price of the moving average option. The same kind of approach 
is used by Grau , but the author improves its numerical efficiency by a different choice of basis 
functions in the regressions used for the conditional expectations estimation. 

Kao and Lyuu [T^] introduce a tree method based on the CRR model to price moving average 
lookback and reset options. Their method can handle only short averaging windows: the numerical 
results that are shown deal at most with 5 discrete observations in the averaging period. Indeed, this 
tree-based approach leads to an algorithm complexity (number of tree nodes) which exponentially 
increases with the number of time steps in the averaging period. Finally, Dai et al. [5] introduce 
a lattice algorithm for pricing Bermudan moving average barrier options. The authors propose a 
finite dimensional PDE model for such options and solve it using a grid method. 

The pricing of moving average options is closely related to high-dimensional optimal stopping 
problems. It is well-known that deterministic techniques such as finite differences or approximating 
trees are made inefficient by the so-called curse of dimensionality. Only Monte Carlo type techniques 
can handle American options in high dimensions. Bouchard and Warin [5] and references therein 
shall give to the interested reader a recent review of this research field. 

More generaly, a related problem is that of optimal stopping of stochastic differential equations 
with delay. With the exception of a few cases where explicit dimension reduction is possible [H [10] , 
there is no numerical method for solving such problems, and the Laguerre approximation approach 
of the present paper is a promising direction for further research. 

The rest of the paper is structured as follows. In Section [2] we introduce the mathematical 
context and provide a general result which links the strong error of approximating one moving 
average process with another to a certain distance between their weighting measures. We then 
introduce an approximation of the weighting measure as a series of Laguerre functions truncated 
at n terms, which leads to (n + l)-dimensional Markovian approximation to the initial infinite 
dimensional problem. The properties of Laguerre functions combined with our strong approximation 
result then enable us to establish a bound on the pricing error introduced by our approach as n 
goes to infinity. In Section[3j our numerical method, based on least squares Monte Carlo algorithm, 
is presented. The final section of this paper reports the results of numerical experiments in the 
Black-Scholes framework which include pricing moving average options with time delay. 

Throughout the paper we assume that the price of the underlying asset S = (St)t>o is a non- 
negative continuous Markov process defined on the probability space (f2, J-",F,P), where P is a 
martingale probability for the financial market and F = (J r t)t<T is the natural filtration of S. 

For the sake of simplicity, we present our results in the framework of a 1-dimensional price 
model but they are directly generalizable to a multi-asset model or to a model with unobservable 
risk factors such as stochastic volatility. 

As usual, we denote by L 2 = L 2 ([0, +oo)) the Lebesgue space of real-valued square-integrable 
functions / on [0, +00) endowed with its norm: 



2 A finite dimensional approximation of moving average options price 
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Strong approximations of moving average processes Consider a general moving average 
process of the forirj^ 



M t 



S t - U fj,(du) 



where \i is a finite possibly signed measure on [0, oo). Throughout the paper, we shall adopt the 
following convention for the values of S on the negative time-axis: 



S t = So, Vi < 0. 



(2) 



We shall use an integrability assumption on the modulus of continuity of the price process: there 
exists a constant C < oo such that 



E 



sup \St — S s 

t,s£[0,T]:\t-s\<h 



(3) 



Fischer and Nappo [9 show that this holds in particular when S is a continuous It 6 process of the 
form 

S t = S + f b s ds + f a s dW s 
Jo Jo 

with 

1 1+7 



E 



sup \b s 

0<t<T 



< oo and E 



SUp \<T S 

0<t<T 



< OO 



for some 7 > 0. 

The following lemma provides a tool for comparing moving averages with different weighting 
measures. 

Lemma 2.1. Let Assumption ^ be satisfied, and let \x and v be finite signed measures on [0, 00) 
with Jordan decompositions \i = /i + — ji" and v = v + — v~ , such that > 0. Define 



M t = / S t - U fi(du), N t = / St-uv(du). 
Jo 



Then 



E 



sup |Af t - JVi| 

0<t<T 



<C|m(R+)-i/(K+)| 



+ C( M +([0,T]) + u-([0,T\) + |/x([0,T]) - v([0,T})\ 



for some constant C < 00 which does not depend on /i and is, where 

F u (t) :=KIM) and F„(t) := /*([0,t]). 



- (4) 



In the literature (see [2] and references therein), moving averages are usually defined via the stochastic integral of 
S. Our definition as an ordinary integral with respect to a weighting measure is closer to the financial specifications. 
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Proof. Step 1. We first assume that fx and v are probability measures. Let 1 and F v 1 be 
generalized inverses of /x and v respectively. Then, 



IE 









r 1 


sup \M t -N t \ 


= E 


sup / \S t 


0<t<T 




0<t<T Jo 






i 

E 


sup | S t 








0<t<T 



-s t 



u)\ du 

du 



L (u) u t-F~ L {u) 

<C I E (\F- l (u)AT-F- l (u)AT\)du 



| J F- 1 (u)AT- J F- 1 (u)AT|du , 



< Ce 



where the last inequality follows from the concavity of e(h). The expression inside the brackets 
is the Wasserstein distance between the measures /i and v truncated at T. Therefore, from the 
Kantorovich-Rubinstein theorem we deduce 



E 



sup \M t — N t \ 

0<t<T 



< Ce 



\FM) - F v (t)\dt 



Step 2. Introduce ft = jul[o,T] an d v — v ^[o.t] + — ^([0, T]))S2t, where 62T is the point 

mass at the point 2T. Then, 



E 



sup \M t 

0<t<T 


-N t \ 


< C|/i(M+) 


- v(R + 


< C|/i(E 


1+) - i/(R+) 






(A+0M + 


0~(R 


+ ))E 


sup 

0<t<T 


Jo 



sup 


f 


0<t<T 


Jo 



ft + (du) + v {du) 
''ft+(R+) + i>-(R + ) ' Jo 



St-u — 



S t - U i>(du) 



ft (du) + i> + (du) 



ft+(R+) + i>-{R+) 



Since ft + (R+) + i> (R+) = ft (R+) + v + (R+) , both measures under the integral sign are probability 
measures, and we can apply Step 1, which gives 



E 



sup |M t -JV t | 

0<t<T 



-C(ft+(I 



< CV(R + ) - kk+)| 



ft+(R+) + D-(R+)J 
C|m(K+) - + C{ft+{R+) + 9~(R + ))e 



\Fn ++i >-{t) - F,- +;>+ {t)\dt 



ft+{R+) + P-(M+) Jo 



\FJt)-F v (t)\dt\, 



because ft coincides with /i and v coincides with v on [0,T]. Using the properties of the function e 
and the definition of ft and i>, we then get Q with a different constant C. □ 



Introducing Laguerre approximation The aim of this paragraph is to provide heuristic argu- 
ments which lead to Laguerre approximation of the moving average. We would like to find a finite- 
dimensional approximation to M, that is, find n processes Y 1 , . . . , Y n such that (5, Y , . . . , Y n ) 
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are jointly Markov, and M t is approximated in some sense to be made precise later by M™ which 
depends deterministically on St, Y 1 , . . . , Y". 

Since M is linear in S, it is natural to require that the approximation also be linear. Therefore, 
we assume that Y — (Y 1 , . . . , Y n ) satisfies the linear SDE 

dY t = -AYdt + l(aS t dt + f3dS t ), (5) 

where A is an n x n matrix, 1 is a rt-dimensional vector with all components equal to 1 and a and 
j3 are constants. Similarly, the approximation is given by a linear combination of the components 
of Y: M n = B Y, where B is a vector of size n and _L denotes the matrix transposition. 
The solution to ^ can be written as 



Y t = e- At 



Y Q + f e- A(t - s h{aS s ds + fidS s ) 
Jo 



or, assuming stationarity, as 

Y t = [ e- A{t - s) l(aS s ds + /3dS s ) and M" = f B ± e - A ^- s h(aS s ds + /3dS s ). 

J —oo J —oo 

Integration by parts then yields: 

M™ = l3B ± lS t + f B ± (a- Ap)e~ A{t ~ s) lS s ds := K n S t + f h n (t-u)S u du, 

J —oo J —oo 

Recalling the structure of the matrix exponential, it follows that the function h n is of the form 

K n k 

M*) = E e ~ Pfct E c " f ' ( 6 ) 

k=l i=0 

where ni + . . . + tik + K = n (K is the number of Jordan blocks of A). Therefore, the problem 
of finding a finite-dimensional approximation for M boils down to finding an approximation of 
the form K n 8 (dt) + h n (t)dt for the measure v. This problem is well known in signal processing, 
where the density h of v is called impulse response function of a system, and h n is called Hankel 
approximation of h. For arbitrary /j, and n, Hankel approximations may be very hard to find, and 
in this paper we shall focus on a subclass for which K = 1, that is, the function h n is of the form 

n-i 

h n (t)=e-v t ^c i t i . (7) 

i=0 

This is known as Laguerre approximation, because for a fixed p, the first n scaled Laguerre functions 
(defined below) form an orthonormal basis of the space of all functions of the form Q endowed 
with the scalar product of L 2 ([0,oo)) which will be denoted by (•, •). See [TB] for a discussion of 
optimality of Laguerre approximations among all approximations of type Q. 

Definition 2.1. Fix a scale parameter p > 0. The scaled Laguerre functions {Lt)k>o are defined 
on [0, +00) by 

Ll{t) = ^2pP k {2pt)e- p \ Vfc>0 (8) 



G 



in which (Pk)k>a is the family of Laguerre polynomials explicitly defined on [0, +00) by 

k 



i=0 



k — i 



(-ty 



V/c > 



(9) 



or recursively by 



Po(t) = 1 

P 1 (t) = l-t (10) 
Pk+i(t) = ((2* + 1 - t)P k (t) - kPk-i(t)),Vk > 1. 

The scaled Laguerre functions (L^)k>o form an orthonormal basis of the Hilbert space L 2 ([0, 00)), 

V(j,fc), (Ul,Li) = 5 j , k . 

Fix now an order n > 1 of truncation of the series. In view of Lemma |2.1| we propose the 
following Laguerre approximation of the moving average process M: 

• Let H(x) = [i([x, +00)). 

• Compute the Laguerre coefficients of the function H: 

A P k = (H,Ll). 

-in-1 



Set H'P{t) = Y2=l ^fc L fe(0 and Kit) = -f t HP(t). In view of Lemma 
can be written as 

n— 1 n— 1 

K = J2<Ll, al= P Al + 2 P ^ K- 

k=0 i=k+l 



A.2 



the function hE 



(11) 



• Approximate the moving average M with 



M^ = {H{0)~H?M)S t + 



+ OO 



h p n (u)St- u du, Vt>0. 



(12) 



Remark 2.1. The approximation proposed in ([12j (and in particular the correction coefficient in 
front of St) is chosen so that the total mass of the weighting measure of the approximate moving 
average M™' p is equal to the total mass of the weighting measure \i of the exact moving average. 
In particular, such an approximation becomes exact for a constant asset price S. 



From definitions ( 11 ) and ( 12 (, it seems natural to introduce n random processes X p '°, . . . , X p ' n 1 
defined by 



+00 



L p k (v)S t - v dv,Vt >0,Vife = 0, 



,n — 1. 



(13) 



They will be called Laguerre processes associated to the process S throughout this paper and are 
related to the moving average approximation by 



= (H(0) - H*(0))St + £ alXf' k , Vt > 0. 



(14) 



fc=0 
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Proposition 2.1. Letn> 1 andp>0. The (n+l) -dimensional process (S, X?' , X?' 1 X?' 11 - 1 ) 
is Markovian. The n Laguerre processes follow the dynamics: 

'dX?>° = L/2pS t -pXf )dt 
dXf' 1 = U2pS t - 2pX? fl - pXf' 1 ) dt 



with initial values 



dxr- 1 = (vzps t - ipEiZo^-pxr- 1 



dt 



Xfr k = S o (-l) fc ^,Vfc>0. 



(15) 



Proof. Immediate, from Equations (13) and ([8| and properties A.l(iJ and A.l(n| 



□ 



Convergence of the Laguerre approximation 

Proposition 2.2. Let Assumption ^ be satisfied, and suppose that the moving average process 
M is of the form 

/>oo 

M t = K S t + / S t - U h(u)du (16) 
Jo 

where Kq is a constant and the function h has compact support, finite variation on M., is constant 

0,T}. Then 

< Cs{rri) 



in the neighborhood of zero and is not a.e. negative on [0, T]. Then the error of approximation (14) 
admits the bound 



E 



sup \M t -M^ p \ 



0<t<T 



Proof. We shall use Lemma 2.1 The measures /i and v are defined by /i(dx) = K S (dx) + h(x)dx 



and v{dx) = (H(0) — HP(Q))b {dx) + hf l (x)dx. Therefore, these measures have the same mass, and 
the first term in estimate Q disapears. In addition, 

\(i([0,T})-v([0,T])\ = \H(T)-H%(T)\, 



which remains bounded by Lemma A. 4 Let us show that ^~([0, T]) is bounded as well. For this it 
is enough to prove that ||(i?^)'||2 is bounded on n. A straightforward computation using Lemma 
IA. 21 shows that 



k-i 



k 



2pH{Q)-2pY J A p l - P Al, 



(17) 



i=0 



where c\ := (h,L p ) are Laguerre coefficients of h. By definition of H% and in (11), this leads to 



k 



2p[H(0)-HP(0)] 
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We have thus: 

ll Will = E Ki 2 ^ 2 E Ki 2 + 2 E \V^p[H(o)-Hm} 2 = o(n-i; 



k<n-l 



k<n-l 



k<n-l 



by Lemma A. 4 and using ^/2p [H(0) — H p (0)] — c p n + pA p l issued from (17). Therefore, there exists 
a constant C < oo, which does not depend on n, such that 



E[ sup \M t - Air\] < Ce ] 

0<t<T V/o h+ (t) dt J ° 



\H{t)-H p (t)\dt\ 



By Cauchy-Schwartz inequality and Lemma |A.4| 

£ \H{t) HP(t)\dt <Vf\\H- flS|| a = Vf (j2 \<\ 2 j = 

from which the result follows using the properties of s and the fact that h is not a.e. negative on 
[0, T] (which means that J Q T h + (t)dt > 0). 

□ 

Approximating option prices The price of the American option whose pay-off depends on the 
moving average M and the price of the underlying is given by 

supE[ < /.(S' r ,Af r )] 
rer 

where T is the set of F-stopping times and </> is the payoff function. It can then be approximated 
by the solution to 



supE[</>(SV,M"> p )] 



(18) 



Proposition 2.3. Let Assumption ^ be satisfied, and suppose that the payoff function <f) is Lip- 
schitz in the second variable and that the moving average process M satisfies the assumptions of 



Proposition 2.2 Then the pricing error admits the bound 



^■pricing \P*i P) 



sup E [<t> (S T ,M T )] - sup E [cf> (S T , M?' p )} 

t£T tET 



< Ce(n 



where C > is a constant independent of n. 
Proof. We have first: 

Vr, E [0 (S T ,M T )} — E[cj) {S T , M" ,p )] + E [0 (S T ,M T ) - (S T ,M?*)} 
=> su P E[0(S T ,Af T )] = sup (E[cj>{S T ,M?< p )} +E[<fr(S T ,M T ) - cj){S T ,M^ p )] 

T T ^ 

< supE [0 (S T , M?< p )] + supE [(f> {S T , M T )-<j> (S T ,M?< P )] 
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In consequence, 

SUpE [(j> (S T , M T )] - SUpE [(f) (S T ,M?< P )] < supE \<t> (S T , M T ) - <j) (S T , M^ p )\ . 



By symmetry and (A2), we get 

supE [4> (S T ,M T )] - supE [</> (S T , M?' p )} 

T T 

and the result follows from Proposition |2.2| 



< supE|M. 

< E 



M n,p\ 



sup \M t -M^ p \ 



0<t<T 



Uniformly- weighted moving average The uniform weighting measure 

fi(dx) = h(x)dx = -t[oj]dx 



□ 



(19) 



satisfies the assumptions of Proposition 2.2 In particular, H (x) = i (5 — x) . From Lemma 



the Laguerre coefficients A d ^ p = (H, L p ) arc related to the Lagucrrc coefficients of h, 



A.2 



At* = (-l)^-i^-T(-lf 



k-1 



k-iJ,p 



(h,L p k ) 2 , 



(20) 



and the coefficients of h can be computed from the values of Laguerre polynomials: 



Sp 



(1 - e- p5 P n (2 P S)) + 2^(-l) fc (l - e- pS P n _ k (2 P 5)) 



k=l 



(21) 



Given the length of the averaging window 5 > and an order n > 1 of approximation (number 
of Laguerre functions), we determine the optimal scale parameter p pt{3, n) as 



Po P t($,n) = argmin||£f - H p \\ 2 = argmm • 



p>0 



p>0 



k=0 



(22) 



Finding an explicit formula to p pt(d, n) does not seem to be possible, but finding a numerical 
solution is easy using the explicit expressions ( 20 ) , ( 21 ) and ( 22 ) . In addition, once p opt is computed 
for a couple (l,n), the scaling property of Laguerre functions ^ gives the value of p opt (5, n) , for 
any S > 0: 



J? op t(<5, n) 



Popt(l,™) 



Table [T] gives the values for p op t(l>"-) for the first 10 values of n computed with an accuracy 
of 10 -3 . Figure [l] (left graph) illustrates the approximation of H by the truncated Laguerre ex- 
pansion Hn° pt ^ nS> for n = 1,3,7 Laguerre basis functions (with 6 = 1). The corresponding error 
\\H — Hn° pt ^ W2 as a function of n is shown in the right graph. The error is less than 5% already 
with n — 3. A simple least squares estimation by a power function gives a behavior in 0(n -1 ' 06 ). 
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n 


1 


2 


3 


4 


5 


6 


7 


8 9 


10 


2>opt(l,n) 


2.149 


4.072 


6.002 


4.234 


5.828 


7.473 


9.155 


10.866 9.153 


10.726 



Table 1: Optimal scaling parameters for approximating H(x) = i (5 — x) + . 



12°o 



1 - 2 1 




0.0 0.5 1.0 1.5 2.0 2.5 3.0 0% 



123456789 10 
n 

Figure 1: Left: Laguerre approximation of the function H(x) = | (S — x) + . Right: L 2 error of the 

approximation. 

Proposition |2.2| also applies when the moving average is delayed by a fixed time lag I > 0: 

1 /■*-' 

X t = - S u du, Vt>6 + 1. (23) 

d Jt-l-6 

In this case, the weighting measure is fj,(dx) = \\[i_i + s\dx and H{x) — \{{8 + 1 — x) + — (I — x) + }. 

Figure [2] shows the approximation of H by Hn° vt ^ for n = 1, 3, 5, 7 Laguerre basis functions with 
5 = 1 and / = 0.5 as well as the L 2 -error made as a function of n (in the same way as above, we 
numerically compute p pt(n) for minimizing the L 2 -error made by Laguerre approximation). In 
comparison to the previous case, it appears that the number of Laguerre basis functions necessary 
to approximate H is greater for an equivalent accuracy of the approximation: the error is less than 
5% from n — 5. This is due to the fact that the density of the weighting measure has two points of 
discontinuity. 
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5 2.0 3.5 3.0 3.5 4.0 



Figure 2: Left: Laguerre approximation of the function H(x) = \ {{8 + I — x) + — (I — x) + }. 

Right: L 2 error of the approximation. 



3 A Monte Carlo-based numerical method 

In this section we present a numerical method for computing the solution to the approximate 



problem (18 1, which is a (n + l)-dimensional optimal stopping time problem. For the sake of 



simplicity, we restrict ourselves to uniformly weighted moving averages. Since the dimension of the 
problem may be high, we use a Monte Carlo technique. Our numerical approach corresponds to 
the one from Longstaff and Schwartz [T5] and the computation of conditional expectations is done 
with a regression based approach. In particular, we shall use the technique of adaptative local basis 
proposed by Bouchard and Warin [5]. 

Forward simulation in discrete time We compute the price of the discrete time version of 



the American option ( 18 ) in which the moving average X has been replaced by its approximation 



M n ' Popt defined in (14) and the exercise is possible on an equidistant discrete time grid with N > 1 
time steps At = 7T = {to = 0, t%, . . . , tjy = T}. We assume that there are exactly Ng > 1 time 
steps within the averaging window of length 5: N$ = -A = 4iV, and that the spot price S can be 



simulated on 7r either exactly or using the Euler scheme. We denote this simulated discrete time 
price by 

extend this definition to [0, T] by 



Sf = Sf , Vte(U,t i+1 ] 



(24) 



and shall also apply convention ^ to S*. We define the discrete time version X* of the moving 
average process X by 



X l = 



u-s 



S?dt = 



1 



^-N s + l 



e 7r. 



(25) 
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Similarly, the discrete time versions of the Laguerre processes are defined by 



Xf^ = t Ll{U-v)S^dt^j2{si-Sl\^ 3 + l)Mct J+mt - V + S {-l) k ^, VUen. 

J-oo j=1 P 

Backward resolution of the optimal stopping time problem The resolution is based on 
the well-known backward American dynamic programming principle. We adopt a Longstaff and 
Schwartz-style approach which consists in estimating the optimal exercise time (or equivalently the 
optimal cashflows generated by the optimal exercise rule) instead of focusing on the computation of 
the option value processes (as for example in Tsitsiklis and Van Roy [H]). Throughout the paper, 
the approach presented below will be called (Lag-LS). The optimal payoffs are evaluated using the 
approximate value of moving average X 71 derived from (14): 



M n,p opt ,« = _ #Popt( ))5£ + a P h ° pt Xl° pt ^,Vt t E 7T. 



(26) 



fe=0 



Denote by (r^)i=ff e if the sequence of discretized optimal exercise times: is the optimal 
exercise time after U £ n. The backward algorithm works as follows: 

1. Initialization: = T 

2. Backward induction for i = N — 1, . . . , N$: 

r? = Ut Ai + r? +1 t CAi with A, = [cf> (SI, M***'*) > E u [<f> (s^M^f"^)] } 

3. Estimation of the option price at time 0: 



v =E 



5^. ,AQ 



»,Popt,ir 



in which: 

E ti [•] = E [.| (S^X^' '*, . . .,X^ pt,n - 1 ' ,r ) 

Estimators of the conditional expectations are constructed with a Monte-Carlo based technique. It 
consists in using M > 1 simulated paths on tt of the (n + l)-dimensional state process: 

^ S n,( m ) jX p opt ,0,n,(m)^ ^ X p opt ,n-l,K,(m)^ ^ m = 1, . . . ,M. 
The corresponding paths of the approximate moving average are denoted by: 

M n,p opt ,7r,(m) j m= l ) ... j M. 

Conditional expectations estimators Ef 1 are then computed by regression on local basis functions 
(see the precise description of the procedure in Bouchard and Warin 5J). We shall denote by 
(b , , . . . , b^-i) the numbers of basis functions used in each direction of the state variable: b for 
, by for X Popt ' 0,7r , bf for X Popt,1 ' ,r , etc. The Monte-Carlo based backward procedure becomes 
thus: 



13 



1. Initialization: t^™ 1 ' = T, m = 1, . . . , M 



2. Backward induction for i = N — 1, 



,iV 5 , m = 1, 



= til .(m) 



7T,(m) 



,M: 



0(s^ +i> Af^-' )]} 



3. Estimation of the option price at time 0: 



•7r,(m) , «-n,p pt,ir,(m) 



Remark 3.1. We will use a numerical improvement to this standard backward induction algorithm, 
which might seem rather natural for practitioners. It consists in evaluating the optimal payoffs using 



the exact value ( 25 ) of the moving average. In particular, the optimal stopping frontier becomes: 

A* = (SI, XI) > E ti [</, (s^ +i ,X- f+i )] } . 

This improved method will be called (Lag-LS*) and unlike (Lag-LS) will exhibit a monotone con- 
vergence as n goes to infinity. 



"Non Markovian" approximation for moving average options Motivated by a reduction of 
dimensionality, the numerical method that is most often used in practice to value moving average 
options consists in computing the conditional expectations in the Longstaff-Schwartz algorithm 
using only the explanatory variables (S, X): namely, the price and the moving average appearing in 
the option payoff. The resulting exercise time is thus suboptimal, but the approximate option price 
is often close to the true price. To assess the improvement offered by our method, we systematically 
compare our approximation to this suboptimal approximate price, also computed using a Longstaff 
and Schwartz approach and referred to as (NM-LS). 

Let (0f)i = N s ,...,N denote the discrete time sequence of the estimated optimal exercise times (0? 
being the optimal exercise time after ij 6 7r). (NM-LS) algorithm works as follows: 



1. Initialization: 6Jj = T 

2. Backward induction for i = N — 1, 



,N S : 



6? = Ul Ai + 9f +1 l CAi with Ai = {0 (SI, XI) > E [cb (s% +i ,X$ f+ J I (SI, XI)] } 



3. Estimation of the option price at time 0: 



US = E 



Similarly to other methods, the conditional expectations are computed with the adaptative local 
basis regression-based technique from [S]. The numbers of basis functions used in each direction 
will be denoted by b s for and b x for X 71 . 
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4 Numerical examples 



For our examples, we use the single-asset Black and Scholes framework. We study standard moving 
average options for different values of the averaging window S as well as moving average options 



with delay ( 23 1 



With the same notations as in Section [3j recall that the dimension of the discrete time version 
of moving average option pricing problem is equal to Ns with a Markovian state: 



We use the standard Longstaff and Schwartz algorithm for such a Bermudan option in dimension 
Ns as the benchmark method. This method will be called (M-LS) and our Monte Carlo regression 
based approach (see more details in [5]) allows to deal with cases up to dimension 8. For applications 
in which Ns is larger, this method becomes computationally unfeasible. 

We provide at the end of this section a numerical comparison between the convergence rate of 
our Laguerre-based approximation and (M-LS) with respect to the state dimension. 



Moving average options: benchmark prices Consider a standard moving average American 
option with value at time 0: 



sup E [e- rT (f>(S T ,X T )] , X T = 




where the asset price S is assumed to follow the risk-neutral Black and Scholes dynamics: 

dS t = S t {rdt + adWt) , S = s 

and W is a standard Brownian motion. We shall consider call options with pay-off <f>(s, x) = (s— x) + . 
Unless specified otherwise, the following parameters are used below: 



Maturity 


T = 


0.2 


Risk free interest 


r = 


5% 


Volatility 


a = 


30% 


Initial spot value 


s = 


100 



and we consider a Bermudan option with exercise possible every day (when T = 0.2, the time 
interval [0, T] is divided into N = 50 time steps). 

Table [2] shows the prices of moving average call options computed by (NM-LS) and (M-LS) 
for various averaging periods S, with M = 10 million of Monte Carlo paths and b s — b x = 2. 
The prices are averages over 5 valuations and the relative standard deviation is given in brackets. 
For reasonable volatility coefficients of the underlying price process and relatively small averaging 
window (5, (NM-LS) seems to provide a very good approximation (from below) to moving aver- 
age options prices. This justifies the approximation made by practitioners and (among others) by 
Broadie and Cao [4]. 
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N 5 


(NM-LS) 


(M-LS) 


2 


1.890 (0.011 %) 


1.890 (0.011 %) 


3 


2.684 (0.011 %) 


2.685 (0.010 %) 


4 


3.183 (0.018 %) 


3.186 (0.012 %) 


5 


3.526 (0.016 %) 


3.531 (0.007 %) 


6 


3.773 (0.016 %) 


3.780 (0.013 %) 


7 


3.955 (0.011 %) 


3.964 (0.215 %) 


8 


4.092 (0.015 %) 


4.103 (0.316 %) 


9 


4.193 (0.016 %) 




10 


4.268 (0.019 %) 





Table 2: Moving average options pricing with (NM-LS) and (M-LS). 



115 



110 - 



80 




— Spot price S 


— Moving average 


Moving average approx. with n = 1 


a Moving average approx. with n = 3 


* Moving average 


approx. with n = 7 



Time step 



Figure 3: Simulated trajectory of the moving average process and its Laguerre approximations. 



Moving average American options: Laguerre approximation Figure [3] shows a simulated 
trajectory of the underlying price S n , its moving average X n with 8 = 0.04 and the corresponding 
Laguerre-based moving average approximation M rl,Popt ' 7r with n — 1,3 and 7 Laguerre basis func- 
tions. Already for n > 3 Laguerre basis functions, Af" ,Popt,7r accurately mimics the exact moving 
average dynamics of X* and this approximation seems to be almost exact when n = 7. 

Table [3] reports the prices of moving average call options computed using the Laguerre-based 
method presented in Section [3] (Lag-LS) and its improved version (Lag-LS*) (with the same pa- 
rameters as above, in particular S = 0.04). The price values are means over 5 valuations, the 
relative standard deviation is given in brackets and we used M — 5 million Monte Carlo paths for 
n = 1, . . . , 3 Laguerre functions and M — 10 million Monte Carlo paths for n — A, ... ,7 Laguerre 
functions, with 6 s = 4 and b£ = l,Vfc > 0. With M = 10 million Monte Carlo paths, b s = 4 and 
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b x = 1, (NM-LS) gives an option value equal to 4.268. 



n 


(Lag-LS*) 


(Lag-LS) 


1 


4.266 (0.020 %) 


4.092 (0.017 %) 


2 


4.273 (0.022 %) 


4.302 (0.019 %) 


3 


4.276 (0.023 %) 


4.182 (0.018 %) 


4 


4.276 (0.022 %) 


4.227 (0.020 %) 


5 


4.277 (0.023 %) 


4.275 (0.020 %) 


6 


4.277 (0.024 %) 


4.287 (0.022 %) 


7 


4.277 (0.024 %) 


4.258 (0.022 %) 



Table 3: Moving average options pricing with (Lag-LS) and (Lag-LS*). 



Remark 4.1. When the averaging window is large, the variance of the Laguerre states (X%' k )k>o 
is small, and at least much smaller than the variance of the price S. In consequence, increasing the 
numbers b x of basis functions in the directions of these states does not have a strong impact on 
the conditional expectation estimators and the resulting option price. On the contrary, the number 
b of basis functions in the direction of the spot price S should be sufficiently large. 

Whereas (Lag-LS) oscillates as n increases (this is due to the non monotone approximation of 
the moving average X by M" ,p °p t ), (Lag-LS*) shows a monotone convergence when increasing n, 
as shown in Figure [4j The limiting value (almost 4.277) is around 0.2% above the value computed 
by the practitioner's approximation (NM-LS). 



4.280 

4.278 

4.276 

a) 4.274 

c 4.272 
o 

"5. 

° 4.270 
4.268 
4.266 
4.264 



■■—(Lag-LS*) 

Benchmark by (NM-LS) 



Figure 4: Convergence of the improved Laguerre-based approximation. 



Figure [5] presents the prices of moving average call options computed by (Lag-LS*) and (NM-LS) 
when varying 5 from to T with the same parameters as above. 7 Laguerre basis functions were 
used with method (Lag-LS) as soon as Ns > 8. For smaller N$, we take n = N$ — 1: n must satisfy 
the condition n < Ng — 1 because otherwise the estimation of the conditional expectation at time 
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£jv a leads to a degenerate linear system. In the limit case when 5 = T, we retrieve the price of the 
Asian option with payoff (jSt — ^ Jq Stdtj . For large averaging periods, the price that we obtain 

with 7 Laguerre functions is about 0.30% above the benchmark value given by the (suboptimal) 
non Markovian approximation. 




Figure 5: Moving average option price as function of the averaging window 5. 



Moving average American options with time delay Consider now a moving average Amer- 
ican option with time delay I > whose value at time is: 



sup E [0 (S T ,X T )} , X T = \ f S u du. 



With the same option characteristics and parameters as above and an averaging period equal to 
S = 0.02 (number of time steps N$ — 5), Figure [6] presents the prices of delayed moving average 
call options computed by (Lag-LS*) and (NM-LS) when varying I from to T — 5. In the limit case 



+ 



when I = T — S we retrieve the price of the Asian option with payoff ( St — Y-i Jq Stdt 

The relative difference between the option values given by (Lag-LS*) and (NM-LS) is significant 
(bigger than 5%) for time lags such that I € [0.04,0.152] (corresponding to 10 < Ni < 38). For 
example, when I = 0.1 (corresponding to Ni — 25 time steps), the relative difference is around 11%. 
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15 20 25 30 

time lag number of time steps 



45 



Figure 6: Price of the moving average option with time delay as function of the lag I. 

Now fix I — 0.08 (Ni = 20). As shown in Figure [7j when the averaging window increases this 
relative difference decreases. But it is still around 5% when N$ = 15. 



7.50 



7.00 



4.00 - 



3.50 



3.00 



(NM-LS) 

(Lag-LS*) with n = 7 
(Lag-LS*) with n = 5 
(Lag-LS*) with n = 3 
(Lag-LS*) with n = 2 




5 10 15 20 25 

number of time steps within the averaging window 



30 



Figure 7: Price of the moving average option with time delay 
as function of the averaging window S. 
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As expected, for pricing delayed moving average options, the non Markovian approximate 
method (NM-LS) seems to be a worse approximation than in the case without time lag: the error 
increases with the time lag but decreases with the length of the the averaging period. For large 
time lags and relatively small averaging period, our improved Laguerre-based method (Lag-LS*) 
gives option prices up to 11% above the benchmark value given by the suboptimal approximation 
(NM-LS) (cf. Ni — 25 on Figure [6]). For a good accuracy of (Lag-LS*), the required number of 
Laguerre functions is however bigger than in the case without time lag as explained at the end of 
Section[2] 

Pricing of moving average Bermudean options: a convergence rate improvement To 

compute the exact price of a moving average Bermudan option, one can either use the classical 
method (M-LS) taking a sufficient number of steps within the averaging window or use the Laguerre- 
based method with sufficient number of state processes. In this last example we illustrate the fact 
that our method (Lag-LS*) converges much faster than classical method (M-LS) with respect to 
the state dimension for pricing the same Bermudan option. 

Let us consider a moving average call option with maturity T = 0.5 and moving window S = 0.1. 
Figure [5] provides a comparison between price values given by (Lag-LS*) for a time step At = 
when varying the number of Laguerre functions from 1 to 7 and by (M-LS) when varying the num- 
ber of time steps within the averaging period from 2 to 8, that is At = . . . , (the state 
dimension varies in both cases from 2 to 8). M — 20 million of Monte Carlo paths were used in 
both cases and b s = 2, b x = 1. 




Figure 8: Convergence of the improved Laguerre-based approximation and the benchmark 

method for pricing a Bermudan option. 
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A Appendix 

Lemma A.l. The Laguerre polynomials (Pk)k>o belong to C°° ([0, +oo)) and: 

(i) Vfc > l,tP' k {t) - kP k {t) + kP k -x{t) = 

(ii) Vfc > 1, f (P k (t) - P k ^(t)) = - E^To Pi(t) 

Proof. (JTJ) can be found for example in Szego [21] and ([h]) is a consequence of (10). □ 

Lemma A. 2. The definite integrals and derivatives of Laguerre functions can be computed using 
the following formulas: 

/ e- s ' 2 P n {s)ds = 2e- t ' 2 P n {t) + 4e~ 4 / 2 ]T(-l) fe P„_ fc (i). (27) 
Jt k=i 

(e-^P n {t))' = -]T e -*/2p fcW _ I e -*/3p n (t). (28) 

Proof. This follows, after some computations, from the contour integral representation of Laguerre 
polynomials: 

1 f e~^ 

Pn(t) = — & 7 : TTds. 

v ' 2ni J (I- s)s n + 1 

□ 
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Lemma A. 3. The Laguerre functions and their integrals admit the following representation in 
terms of Bessel functions: 

oo ^2 

e-^ 2 P n (x)=J2 A k(^) Jfc(v^i). ^ = 4n + 2 (29) 



k=0 



I n (x) := / e~ x l 2 P n (x')dx' = 2^A fc (-) J fe+ i(^) 

;. n ^ 



(30) 
(31) 



where A = 1, Aj = 0, A 2 — \ and other A r s satisfy the equation (to + 2)A m+2 = (m + l)A m — 
^A m -i. The series converge uniformly in x on any compact interval. 

Proof. The first formula is from [7 . The other two follow readily using the integration formula for 
Bessel functions: 

r-l 

X u+ J v {ax)dx = a~ 1 J u+ i(a). 

□ 

Lemma A. 4. Let /i be a finite signed measure on [0,oo), with bounded support which does not 
contain zero. Let {c n } denote the Laguerre coefficients of the function h(x) := /j>([x,oo)) and {A n } 
denote the Laguerre coefficients of the function H(x) := h(t)dt. Then 

Cn = 0(n- 3/4 ) and A n = 0(n- 5/4 ). 
In addition, for x > fixed, e- x / 2 P n {x) = (Din- 1 / 4 ). 

Proof. This result follows from Lemma |A.3[ using the asymptotic expansion for Bessel functions 

J n {x) = { l -vxy^ cos (x - \n + 0(x" 3 / 2 ), 
which holds uniformly ,7j on bounded domains outside a neighborhood of zero. □ 



23 



